subroutine ENERGY(mass,V,psi,E)
!calculates the energy

parameter (nx = 10001, NSTATES = 8)
double complex :: psi(NX,NSTATES), CKIN(NX), CA(NX), CB(NX), CC(NX)
real*8 mass(0:nx+1), SUMV, SUMK
real*8 E(NSTATES), V(NX)
integer ix, is
double precision :: dangle, Ry, a0, am, am1, am2, ec, pi, radius
double complex :: dt, img

common / constants / dangle, dt, img, Ry, a0, am, am1, am2, ec, pi, radius

dx = radius * dangle
!write(*,*) dx
!stop

do is = 1,NSTATES

	SUMV = 0.D0
	do ix = 1, nx
 		SUMV = SUMV + DCONJG(psi(ix,is))*V(ix)*psi(ix,is)*dx
	enddo

	SUMK = 0.D0
	do ix = 1, nx
		if(ix.eq.nx)then
			CA(ix) = 0.D0
		else
			CA(ix) = -psi(ix+1,is)*(1.D0/mass(ix+1) + 1.D0/mass(ix))/2.D0/dx**2
		endif

 		CB(ix)=psi(ix,is)*(1.D0/mass(ix+1) + 2.D0/mass(ix) + 1.D0/mass(ix-1))/2.D0/dx**2

		if(ix.eq.1)then
			CC(ix) = 0.D0
		else
			CC(ix) = -psi(ix-1,is)*(1.D0/mass(ix) + 1.D0/mass(ix-1))/2.D0/dx**2
		endif

		CKIN(ix) = CA(ix) + CB(ix) + CC(ix)

		SUMK = SUMK + DCONJG(psi(ix,is))*CKIN(ix)*dx
	enddo

	E(is) = (SUMV + SUMK)*Ry
enddo

return
end
